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Abstract 

Semi-Markov processes (SMPs) provide a rich framework for many real-world problems. 
However, due to difficulty implementing practical solutions they are rarely used with their 
full capability. The theory of SMPs is quite mature but was mainly developed at a time when 
computational resources were not widely available. With the exception of some of the sim- 
plest cases, solutions to SMPs are inherently numerical, and SMPs have been underutilized 
by practitioners because of difficulty implementing the theory in applications. This paper 
demonstrates the theory and computational methods needed to implement SMP models in 
practical settings. Methods are illustrated with an application modeling the movement of 
coronary patients in a hospital. Our aim is to allow practitioners to use richer SMP models 
without being burdened with the rigorous mathematical theory. 
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1 Introduction 



Semi-Markov processes (SMPs) provide a rich framework for many real-world problems. SMPs 
include Markov processes, Markov chains, renewal processes, Markov renewal processes, Poisson 
processes, birth and death processes, and M/G/l queues to name a few. They have been ap- 



plied in the areas of survival analysis (Butler and Huzurbazar 1997), reliability (Limnios and 



Opri§an 


2001 


), DNA analysis ( 


Barbu and Liminios 


nance ( 


Janssen and Manca 


2007 


), and many others. 



2008), queuing theory (Kleinrock 1975), fi- 



A key reason for this is the perceived complexity of 
models such as continuous-time 



as one would expect, given their generality 
solving semi-Markov models. Many practitioners prefer "simpler 
Markov chains, but (as we demonstrate in this paper) solving semi-Markov models for quantities of 
interest is in fact straightforward. Within the realm of finite-state processes, SMPs are extremely 
flexible, and we believe they should be a primary statistical tool for modeling stochastic processes. 
The concept of semi-Markov processes is generally agreed to have been independently intro- 



duced by Levy (1954), Takacs (1954), and Smith (1955). The theory was formalized in Pyke 
(1961a) and Pyke (1961b). Further developments came from Takacs (1959), Pyke and Schaufele 



(1964), Pyke and Schaufele (19661), and Cinlar (1969). Since that time there have been many 



contributions to and applications of SMPs; for typical applications, however, we feel that the key 
problems in SMPs were solved many years ago, and the relevant theory is treated comprehensively 
in the cited papers by Pyke. Pyke's work provided general solutions in terms of Laplace transforms, 
however, the actual quantities of interest could not be obtained without numerical computation to 
recover probability distributions. The necessary numerical algorithms and computational power 
needed to solve SMPs are readily available today (as we will show), but were not within reach for 
most researchers in the early 1960s when the major results for SMPs were published. Therefore, 
the bridge from theoretical solutions to practical implementation was not immediately established. 
Statistical flowgraph models (SFGMs), as detailed in Huzurbazar (2005), were a beginning in this 
endeavor. SFGMs have primarily focused on finding the first passage distribution from a beginning 
state to some absorbing final state, using methods developed in Mason (1953). This paper more 
fully explores connecting the theory of SMPs with computational methods. In particular, solutions 
for several key quantities found in Pyke (1961b) are implemented using the R software program 
without the need for any additional packages or add-ons. However, these calculations could be 
performed in any program with numeric integration and matrix manipulation capabilities, e.g. 
MATLAB, Mathematica, and FORTRAN. 

So far we have used the term "solve" ambiguously; more precisely, it is the process whereby we 
calculate certain quantities from an SMP. Some of the quantities desired from a time-homogeneous 
finite-state SMP, given the process was started in state i at time t = 0, are: 



1. the probability the process is in state j (as a function of time), 

2. the first passage distribution of the time to reach state j, 

3. the probability of reaching a state j, k number of times (as a function of time), 

4. the probability of reaching a state j, k or fewer times (as a function of time), 

5. the expected number of times the process has been in state j (as a function of time), 
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6. the long run probability that the process will be in state j. 

These are significant quantities that modelers find difficult to compute for even simpler Markov 
processes, but the theory and computational techniques currently exist to solve them in general 
terms for SMPs. With today's computational resources, processes that have many states with 
any smooth transition distributions can be solved in a matter of seconds. In the remainder of the 
paper, we detail how to find these quantities. 

In our experience, the primary focus in stochastic processes has been finding asymptotic prob- 
abilities. However, there is a considerable need for finding time-dependent state probabilities of a 
process, before it is asymptotically "well behaved". For example, if we were developing a business 
model it would be very reasonable for investors to want to know the probability of bankruptcy 
at 1 month, or 2 years into the future. These probabilities, as time goes to infinity, are of little 
consequence. 

As another example, the reliability model in Figure [I] has a single absorbing state (state 3), so 
the asymptotic probabilities are trivial: 1 for the absorbing state, and for the other (transient) 
states. In the context of this reliability example it would be very valuable to find the six numbered 
quantities listed above, i.e., to know how the process behaves before reaching state 3. If the system 
happened to be an automobile, then it would be worth knowing the probability that a vehicle is 
in the "Working" or "Under Repair" states, as well as calculating the first passage distribution 
to the "Unrepairable" state. Other items such as the probability of k repairs before time t or 
the expected number of repairs up to time t could be used for developing efficient manufacturer 
warranties. 




Working Under Repair Unrepairable 



Figure 1: A reliability semi-Markov process. 



2 Background 

A real deterrent to using SMPs is the effort required to decipher the notation. In this section we 
define the notation used in this paper as a reference for the reader, and list assumed properties 
of the models we deal with. We have attempted to conform to the SMP literature as closely 
as possible. We also introduce Laplace transforms, give some of the properties that make them 
useful, and show how to numerically invert them into probability density functions (PDFs) or 
cumulative distribution functions (CDFs). We avoid going into measure theoretic details to make 
the paper as accessible as possible, so we assume absolutely continuous distributions. However, 



the SMP results here have been proven using general probability measures; see Pyke (1961a). 
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2.1 Notation and assumptions 

We assume a finite-state SMP with n states running in continuous time. Roughly, this means 
that transitions between states constitute a Markov chain with transition probability matrix [p^] ; 
however, waiting times (times between transitions) may have arbitrary distributions depending 
on both the origin and destination states for the transition. This may be clearer if we contrast an 
SMP with a discrete-state Markov process: once the current state of a Markov process is known, 
no additional information about its future is provided by its past history; this implies that waiting 
time distributions must be memoryless (geometric in discrete time, or exponential in continuous 
time), and independent of the destination state. In a semi-Markov process, since waiting time 
distributions may be arbitrary, knowledge of how long the process has waited in its current state 
can provide information on how much longer it will wait, and what the destination state will be. 

We also make the following assumptions: the process is time-homogeneous, i.e., the waiting 
time distribution for a transition i — > j does not explicitly depend on how long the process has 
been running; there are no instantaneous transitions from a state, which implies that the number 
of transitions in any finite time period is finite; and there are no transitions from a state to itself. 
The last assumption can be made without loss of generality by transforming to an equivalent SMP 



without self-transitions (Pyke 1961a, Lemma 3.2). 

Important classes of states in SMPs are absorbing states, which, once entered, are never left; 
transient states, which with probability one are eventually left and never returned to; and recur- 
rent states, which with probability one are eventually entered, starting at any state in the process. 
Many applications in survival and reliability analysis have one or more absorbing states, repre- 
senting failure or death, with all other states transient. Other application areas, such as queuing 
and inventory processes, typically have all states recurrent. 

The random variable X represents occupancy time in a state, as distinguished from T, the 
time since the process was started or "calendar time" . Most frequently we are interested in finding 
functions of the calendar time t, and unless explicitly stated, when we use the word "time", we 
are referring to calendar time. Next, let Z(t) denote the state of the SMP at time t. Z(t) can take 
on any positive integer value up to and including the number of defined states in the process. 

Under the assumptions given above, there is a one-to-one correspondence between the SMP 
Z{t) and an equivalent Markov renewal process (MRP) N(t) = (Ni(t), . . . , N n (t)), where each 



N.(t) counts the number of transitions into state j that have occurred up to time t (Pyke 1961b). 
It turns out that it is more natural to think of many of the quantities associated with SMPs in 
terms of the counting processes N.(t). For an introductory discussion of MRPs in relation to 



SMPs, see Ross (1992). 



The table below shows the items necessary to define a SMP: 



n 


the number of states in the SMP 




the CDF of the waiting time distribution 


from state i to state j 




the probability that the next state in the 


process is j, given the 




process entered state % 





Constraints on p.. are Y^j=\Pa = 1 ( or ec L ua l to if i is an absorbing state), and all p tj > 0. 
A boldface symbol without indices refers to a matrix valued function (e.g., F(t) is a matrix with 
elements F. . (t) in the i th row and } th column). The Laplace Transform (LT) of a function is denoted 
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with a tilde, so, e.g., F.(s) is the LT of F t .(t). Similarly F(s) is the matrix valued function where 
the LTs are of the individual elements. We will sometimes drop the (x), (t), or (s) on functions 
to reduce notational clutter. For example, one can read F as F(s). It is worth mentioning that 
F(x) is in general not symmetric, so F i .(x) 7^ F .(x). 

Now that we have detailed what is needed to define a SMP, we formally define the quantities and 
functions we want to obtain from a SMP. These are listed in the same order as in the Introduction 
(p. [3]) with a few additions; they are all conditioned on a given starting state i: 



PS) 


LllC UlUUcLUlllLV Lilt- UlUL/tioo lo 111 oLcLLc J cLL Lllllc L 




p(N j (t)>a\z{o)=i) 

the distribution of the time of first passage to state j 


*>y(M) 


P(N i (t) = k\Z(0) = i) 

the probability of reaching a state j, k number of times by time t 




P{N s (t) < k\Z(0) = i) 

the probability of reaching a state j, k or fewer times by time t 




J5[JV,(t)|Z(0) = i] 

the expected number of times the process has been in state j at time t 


71.. 


\im t ^ oo P(Z(t)=j\Z(0)=t) 

the long run probability that the process will be in state j 


Z(t) 


An integer value denoting the state of the process at time t 




The number of times the process transitioned to state j since t = 



For notational convenience, when the derivatives exist we also define: 



/«(*) = 


—F..(x) 
dx ljV ' 


(1) 




^G.m 

dt vKJ 


(2) 




Pij J~ ij 


(3) 


-{! 


if i j 
if % = j 


(4) 



= ( 5 ) 

I is defined to be the n x n identity matrix, and 1 to be the n x n matrix where every element 
is 1. The matrix operation "o" such that [AoB] = A^B^ is the elementwise product of two 
matrices (the Hadamard product). 

2.2 Laplace transform inversion 

In the statistical community there has been considerable avoidance of transforms, even though 
transforms (plus the computations involved in moving between the time and transform domains) 
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are what actually simplifies the problem of solving SMPs by making it similar to solving a system 
of linear equations. 



A prevalent opinion among researchers is exemplified by this quote from |Janssen and Manca 
. . it is then necessary to invert Laplace transforms to obtain the solution of the 



(2006): 



equation and unfortunately, it is well known that this inverse transform is not stable numerically 
. . ." (p. 92); they cite Bellman et al. (1966) and Csenki (1994) as support. Indeed, it has 



been well-known at least since the book by Bellman et al. that the general problem of Laplace 
transform inversion is ill-posed, and the difficulty of inverting transforms of general functions has 



been much discussed since then (Davies and Martin; 1979 Duffy 1993). 



What the conventional wisdom misses is that in analyzing SMPs we are interested only in 
inverting transforms of probability densities and distribution functions of non-negative random 
variables, functions with very special properties, and it has long been known that this restriction 
significantly simplifies the problem, particularly in the case of absolutely continuous distributions 
(Dubner and Abate; 1968; Abate and Whitt 1992 Strawderman; 2004). Our own experience 



confirms that methods based on Abate and Whitt ( 1992 , 1995 ) are strikingly accurate for numerical 
inversion of transforms specified either in closed form or numerically. These methods are also 
simple to implement and require minimal tuning, in contrast to algorithms such as the Weeks 
method, based on orthogonal expansion in Laguerre polynomials (Csenki; 1994, Sect. 10.2). 
Moreover, LTs need not be known in closed form if we are able to compute them by numerical 
integration. 

The Laplace Transform of a generic function </>(£) is defined as: 



<p(s) 



e y>(£) d£, where s is a complex number. 



(6) 



A useful identity for numerically calculating LTs is Euler's formula (not to be confused with the 
EULER method in this paper): e ld = cos(6 l ) + ism(9), where i = yf—1. If we let s = x + iy, 
Equation ^ then becomes: 



<p(s) 



poo poo 

Jo Jo 



(7) 



which enables us to calculate the LT without any numeric complex integration. The inverse 
Laplace Transform can be defined by the Bromwich contour integral, where the complex contour 



can be any vertical line s = a such that <p{s) has no singularities on or to the right of it, see Abate 



and Whitt (1995). This integral is: 



1 

2m 



a+ioo 



ip(s) ds 



(8) 



Again we emphasize that we are avoiding some measure theoretic details, so although for any CDF 
F, ftfix) may not exist, the LT f ii (s) = J °° e~ sx dF(x) is always defined and uniquely represents 



that random variable (Feller 1971| Chapter XIII). 

A fundamental property of the LT, and the basis of its usefulness for solving stochastic pro- 
cesses, is the fact that it maps convolution to multiplication. Thus if the random variables X and 
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Y have densities fx and fy, the LT of the density of their sum is fxfy- This allows complex inte- 
gral equations to be solved algebraically in the transform domain. Since the density f(x) = F\x), 
another well-known property of LTs that is used extensively in this paper is 



F(s) 



(9) 



LTs have very useful properties, but these come at a price. Many parametric distributions do 
not have closed-form LTs, and even if they did, after manipulation in the LT domain they most 
likely would not have closed-form inverses. With today's computational power and algorithms, 
however, this is no longer the significant obstacle it once was. 

There are various ways to numerically invert LTs. The methods we apply in this paper work 
well in terms of speed and accuracy for most distributions; this is an active area of research, and 
we expect to see improvements in the future. In this paper we use the EULER method as shown in 
Abate and Whitt (1995). EULER is specially suited for inverting LTs of probability distributions 
that are smooth (i.e., the CDF has at least two continuous derivatives). Therefore to accurately 
invert LTs of non-smooth distributions other techniques should be considered. Other possibilities 



include using a modification to the EULER method such as in Sakurai (2004), or using a discrete 



approximation with the use of the inverse discrete Fourier transform. 



In simple terms, the EULER method approximates the exact LT inversion integral (Doetsch 



1974) by the trapezoidal method; the result is a Fourier cosine series approximation, and by a 



careful choice of the step size the cosine terms become ±1. The series is then summed using Euler 
summation, which gives the method its name. For a function ip(t) and its LT <f(t), the result is 



<p(t) 



N 

E 

j=Q 



wj Re 



_ / A jiti 

n 



2t 



+ 



(10) 



where Re[-} gives the real portion of a function and Wj is a weight associated with each term (for 
Euler summation the weights are binomial coefficients). N, A, and the u^'s control the accuracy 
of the approximation; see Abate and Whitt ((1992) for details. 



3 Semi-Markov Process Theory 

Now that we have introduced the terminology and methods of computation we move to the main 
focus of the paper. Figure [T] shows a simple example of a semi-Markov process that might be 
encountered in a reliability application. States 1 and 2 are transient, and state 3 is absorbing. 
Here we might be interested in quantities such as the number of times the process is in state 1 in 
a 5 year period, the distribution of the time until the process reaches state 3, or the average time 
spent in state 1 from to 3 years. In this section we show how these questions and others can be 
answered, by finding the LTs of the quantities P i:j (t), G^it), v i .(k;t), V i .(k;t), and M i;j {t). 

3.1 Time-dependent state probabilities 

The time-dependent state probabilities P l3 {t) are invaluable pieces of information from a multi- 
state model. For example, in a reliability application, time-dependent state probabilities could 
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provide probabilities of being in an "Unrepairable" or "Working" state at a given moment in time. 
At the beginning of a process these probabilities might be significantly different from the long run 



probabilities. The formula to find these probabilities is found in Equation (4.1) in Pyke (1961b), 
which is 



P(s) = -(I-q(s)j (I-hOO). (11) 

Given the LTs of these probabilities, the probabilities themselves can be found by inversion back 
to the time domain. The actual computation will be demonstrated in the application section. 

Another quantity of interest that is related to P i (t) is the expected amount of time the process 
will be in a state during a specific period of time. For example in the reliability example in Figure 
[l] it would be of interest to determine the expected amount of time the system is in "In repair" 
during some time interval. The expectation of time in state j in the interval [0,ti) is J* 1 P i:j {t) dt. 
The LT of this integral can easily be expressed as: 



?0 



q( S ) I-hOO). (12) 



Therefore when Equation (12) is inverted at time t\ we obtain j^ 1 P^(t) dt, which is the expected 



amount of time spent in state j from time t — to t — t\ (assuming the process began in state %). 



3.2 First passage distributions 



Other important items available from an SMP are first passage distributions, G ij {t). In the 
reliability example this would give us, for example, the distribution of the time to reach the 
"Unrepairable" state. G(t) is the matrix valued function of the CDFs of the first passage times 
from state i to state j, and g(t) is the corresponding matrix of PDFs. G^(oo) could be less than 
one if there is positive probability that the process may not ever reach state j from state i. These 
distributions are useful in determining how likely transitions are, and how long they may take. 
They are also useful for calculating other quantities of interest. From Equation (4.6) in Pyke 



(1961b) we obtain 



(s) = q(s)(l-q(s)y 1 Io(l-q 



-i 



and we know, 



G(s) 



(13) 



(14) 



At first glance Equation (13) may look intimidating, but for any matrix A, the operation Io A 



just zeros out all the non-diagonal elements of A. Then the inverse of a diagonal matrix is just 
the matrix of reciprocals of the diagonal elements (if they are all nonzero). So once (I — q(s)) 



is calculated, the rest of Equation ( 13 ) is simple matrix manipulation. 



For processes with an absorbing state, such as Figure [TJ we may be interested in finding the 
first passage hazard function. In its basic form the hazard function is defined as: 



X(t) 



f(t) 



F(t) 



(15) 



S 



For the SMP in Figure [TJ the hazard function of the first passage can be calculated by finding 
G 13 (t), g 13 (t), and using Equation (15). 

The first passage distributions can be used to find other quantities, as we see in the next 
section. 



3.3 Markov renewal process probabilities 

The probability v tj (k; t) that the process has reached a particular state a certain number of times, 
is used when the Markov renewal process viewpoint is taken. It provides information about the 
probability that state i has been visited exactly k times by the process up to time t. Again 
considering the reliability example in Figure [TJ this would tell us how probable it is that 10 repairs 
would be needed by a particular time. The cumulative version, how probable it is that 10 or fewer 
repairs would be needed by a particular time, is given by V i .(k; t). The formulas are simpler if not 



presented in matrix form. The first of these is Equation (5.4) of Pyke (1961b) 



vJk;s) 
















9iM 





9ij 



if k e {1,2,3,...} 
if k = 



(16) 



Since these are probabilities of mutually exclusive events, we can add them to find the cumu- 



lative probability. We find that ^2 

%(k;s) = - 

s 

In matrix form this equation is: 

v(M = - 



n=0 ^ij 



[n; s) is a telescoping sum and obtain: 



Ms) 



for k 6 {0,1,2,...}. 



(17) 



S o 



l(lo[g( s )] fc ) ), for fee {0,1,2,...}. 



Presenting these probabilities in terms of first passage times has an intuitive interpretation: 
V i .(k; t) is just 1 minus the CDF of the first passage from state i to state j convolved with the first 
passage CDF from state j to j, k times. In other words the probability that we must repair the 
system k times or less is 1 minus the CDF of reaching the repair state for the first time convolved 
with reaching it k more times. 

The quantities discussed in this section can be very useful for the planning and allocation of 
resources. Viewing v.. (k; t) as the probability mass function of the random variable K, the next 
section discusses the expectation of K. 



3.4 Expected values of the Markov renewal process 

M i .(t), the expected number of visits to state j at time t, is just the expectation of K or 
Ylk^LokViiik-jt). This also can be presented in terms of first passage distributions, but is more 



elegant as presented in Equation (5.11) from Pyke (1961b): 



Mis) 



-i 



(19) 
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M(t) is what defines a Markov renewal process (MRP) and Equation (19) is a theoretical link 
between SMPs and MRPs. Roughly it says that given a SMP, using Equation (19) we can find a 
corresponding MRP and vice-versa. 



3.5 Asymptotic state probabilities 



The final quantities we discuss are the asymptotic state probabilities 7r y . In some instances 7r y 
may not be well defined; this can occur if the parameterization of a transition distribution has an 
infinite first moment (e.g., the Levy distribution). 



According to Equation (7.5) of Pyke (1961b), we have for continuous time distributions: 



7T = lim P..(t) 



G, 



oo 



JIT Xh n ( X ) ^ x 



t— 'TOO 



Jo°° x 9jj ( x ) dx 



(20) 



In words, Equation (20) says that the asymptotic state probability is equal to the probability of 



reaching state j from i, multiplied by the average time in state j, and divided by the average time 
to return to state j, if just arriving in state j. If the expectations of the transition distributions 
are known, then J °° xh u (x) dx is simple to calculate. Usually G tj (oo) can be approximated by 
finding G i] {a), for a a relatively large number. The main difficultly is finding the mean of the 
first passage distribution, J °° xg^{x) dx. If G^ipo) = 1, one way to find the expected transition 
time of g Sj is to use the moment property of the LT, by calculating — j- s g n (s)\ a _ n . If G..(oo) < 1, 
then f° 
G 



xg.. (x 



J 33 

jr.. = Lr .(oo), and if j is a transient state, n.. 

Having introduced the primary equations for SMPs, we now proceed with an application. The 
difficulties in applying the equations from Pyke ( 1961b[ ) are the computations involved in passing 
to and from the LT domain, and the calculation of (I — q(s)) . 



dx 



oo, and Equation (|20|) no longer holds; however, if j is an absorbing state 

0. 



4 Application: Movement of Coronary Patients 



This section focuses on applying the theory from the previous sections. To illustrate, we use 



a semi-Markov process from Kao (1974) that models the movement of coronary patients in a 



hospital, specifically the myocardial infarction (heart attack) positive patients. This general type 
of model for patient movement is of interest to hospital administrators for planning and allocation 
of physical and personnel resources. However, it is quite analogous to other logistic problems such 



as the flow of aircraft through the depot maintenance process as discussed in Huzurbazar and 



Peters (2007). 



Inferential analysis of the movement of coronary patients, such as the selection of transition 



distributions and parameter estimation from the data, can be found in Kao (1974). For this data 



maximum likelihood estimation was used, but Bayesian inference and other methods could also 
be implemented. We do not discuss the inference in this example in order to focus on solving the 
SMP after parameter estimates are obtained. 

Coronary patients are assumed to be in one of nine possible states: coronary care unit (CCU), 
post-coronary care unit (PCCU), intensive-care unit (ICU), medical unit (MED), surgery (SURG), 
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ambulatory care (AMB), extended care facility (ECF), HOME, and DIED. The last three states 
are considered to be absorbing. 

For the matrices below, the rows represent the instance of starting the process in a particular 
state and the columns present a quantity of a state given that row. The probability transition 
matrix of the process, p, is: 



ecu 

PCCU 

ICU 

MED 

SURG 

AMB 

ECF 

HOME 

DIED 



ecu 

/ 0.0000 
0.0192 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

\ 0.0000 



PCCU 
0.7447 
0.0000 
0.5833 
0.0135 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 



ICU 
0.0084 
0.0137 
0.0000 
0.0405 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 



MED 

0.1339 

0.0247 

0.1667 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



SURG 
0.0042 
0.0027 
0.0833 
0.0135 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 



AMB 

0.0063 

0.0027 

0.0000 

0.0270 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



ECF 
0.0000 
0.0577 
0.0000 
0.0811 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 



HOME 
0.0063 
0.8298 
0.0000 
0.7028 
1.0000 
1.0000 
0.0000 
0.0000 
0.0000 



DIED 
0.0962 \ 
0.0495 
0.1667 
0.1216 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 J 



Waiting time distributions for the process are assumed to be Weibull, as defined below. If X 
has a Weibull(7,0) the PDF of X is: 



/*(*) 



l T i-^ P -(xye) 
6 



The transition distributions as given in Kao (1974) are: 

• f x is Weibull(4.738025, 4566277818.13) 

• / a is Weibull(2.207438, 14541.6089) 

• / 3 is Weibull(0.766338, 16.6991) 

• f A is Weibull(2.303331, 1017649.5158) 

• U is Weibull(l. 623492, 4707.3132) 
The transition matrix f(£) is: 





ecu 


PCCU 


ICU 


MED 


SURG 


AMB 


ECF 


HOME 


DIED 
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fx 


fx 
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fx 


fx 
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fx 


fx 
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MED 
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h 
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SURG 
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AMB 
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ECF 





























HOME 





























DIED 


V o 


























(21) 
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With n = 9, p, and f(x) we have defined an SMP. For a given time t we want to find some 
quantity, for example, the value of the first passage CDF from state 1 to state 3 at time t = 10 
i.e., G 13 (10). First we compute f(s) at the points needed by EULER, by elementwise numeric 



integration since the distributions do not have closed-form LTs. Then using Equations (13) and 



(14), we calculate G(s). Finally, the EULER algorithm gives an estimate of G(10), from which 



we take the element on the first row and third column. The computation time to calculate the 



matrix G(10) is less than a second. Similar methods are demonstrated in more detail in Warr and 
uzurbazar (2010). 

For processes with only a few states, we can afford to compute (I — q) -1 many times with 
only a small computational penalty; for processes with many states it would be more efficient to 
compute the inverse symbolically once, using a computer algebra package such as Mathematica 
or Maple, at the beginning of the calculations. 

The time consuming calculation for this SMP is finding f(s). To reduce the computation 
time we only find this once for each time point t and make all of the calculations from it before 
discarding the values. Programming this way makes the code more complex, but speeds up the 
computations. 

For this SMP we calculated the transition probability matrix P(£) at 120 points, spaced at 
12 hours apart. If we assume the patient began in state CCU we find the following probabilities 
in Figure [5J At any given time (a vertical line), the proportion of color of the line represents 
the probability a patient is in the corresponding state. Clearly as time progresses the patient is 
absorbed into ECF, HOME or DIED. 



03 
-Q 
O 




30 

Time in Days 



Figure 2: The probabilities of being in a state at a particular time, given the patient began in the 
CCU state. 



We also obtain first passage distributions from these calculations; from Equations (13) and 
(14) we obtain both g(t) and G(t). Using these we can find a conditional first passage hazard 
using the equation 
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KM 



GJoc) - GJt) 



(22) 



In Equation (22) the usual 1 in the denominator has been replaced with GL.(oo). This makes 



the assumption that the process eventually ends in state j. Therefore Equation (22) gives a 



conditional hazard function, given the process reaches state j. Figure [3] shows the conditional 
hazard functions for the three absorbing states (again assuming the process began in state CCU). 



Hazard for Passage to ECF 



Hazard for Passage to Home 



Hazard for Passage to Dead 




10 20 30 40 
Time in Days 



20 30 40 
Time in Days 



20 30 40 
Time in Days 



Figure 3: The conditional hazard functions for the three absorbing states if the process began in 
state CCU. 



Next we look at some of the functions of the the Markov renewal process associated with this 
SMP. The probability of revisiting a state in this SMP is small, so we only consider a few values of 
k for v(fc; t). We focus on v(fc; t) since V(/c; t) is easily calculated from v(k; t). For the interesting 
cases, the probabilities of not visiting a state before t = 60, v(0; 60) = P(iV.(60) = 0|Z(0) = i), 
are given in the following matrix, where "60" is 60 days. 



v(0;60) = 





CCU 


PCCU 


ICU 


MED 


SURG 


AMB 


ECF 


HOME 


DIED 


CCU 


/ 0.9854 


0.2454 


0.9749 


0.8420 


0.9894 


0.9872 


0.9426 


0.2191 


0.8405 \ 


PCCU 


0.9806 


0.9766 


0.9848 


0.9698 


0.9955 


0.9963 


0.9385 


0.1214 


0.9412 


ICU 


0.9886 


0.4105 


0.9844 


0.8158 


0.9112 


0.9933 


0.9506 


0.2795 


0.7772 


MED 


0.9993 


0.9627 


0.9593 


0.9922 


0.9829 


0.9727 


0.9163 


0.2186 


0.8687 


SURG 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


0.0000 


1.0000 


AMB 


\ 1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


0.0000 


1.0000 / 



The higher the probability in the above matrix, the less likely that the state will be visited. 
The next matrix shows P(N (60) = 1|Z(0) = i) or v(l; 60) for the non-absorbing states. 



v(l;60) = 





CCU 


PCCU 


ICU 


MED 


SURG 


AMB 


ECF 


HOME 


DIED 


CCU 


/ 0.0144 


0.7370 


0.0248 


0.1568 


0.0106 


0.0128 


0.0574 


0.7809 


0.1595 \ 


PCCU 


0.0191 


0.0229 


0.0150 


0.0300 


0.0045 


0.0037 


0.0615 


0.8786 


0.0588 


ICU 


0.0112 


0.5760 


0.0153 


0.1828 


0.0888 


0.0067 


0.0494 


0.7205 


0.2228 


MED 


0.0007 


0.0365 


0.0401 


0.0077 


0.0171 


0.0273 


0.0837 


0.7814 


0.1313 


SURG 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 


AMB 


\ 0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


1.0000 


0.0000 J 
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This matrix shows which states are likely to be visited once, the higher the probability the 
more likely to be visited. This matrix clarifies that even though the process began in a state i, 
AT(i) = until it is visited again. Therefore, AT(t) is counting the number of times the process 
transitioned to state i after t = 0, not the number of times the process has been in state i. 

The next matrix shows v(2; 60) for all the relevant cases. Clearly, it is quite unlikely that any 
of the states will have 2 visits from the same patient. 



v(2;60) 



CCU PCCU 

CCU / 0.00020 0.01722 

PCCU 0.00028 0.00051 

ICU 0.00013 0.01327 

MED \ 0.00001 0.00074 

Next we show the expected number of visits, Mij(t), 
the values if starting in a non-absorbing state. 



ICU MED 

0.00038 0.00120 \ 

0.00023 0.00022 

0.00022 0.00138 

0.00061 0.00005 / 

at 60 days, or M(60). We display only 





CCU 


PCCU 


ICU 


MED 


SURG 


AMB 


ECF 


HOME 


DIED 


CCU 


/0.015 


0.773 


0.026 


0.159 


0.011 


0.013 


0.057 


0.781 


0.159 \ 


PCCU 


0.020 


0.024 


0.015 


0.030 


0.005 


0.004 


0.061 


0.879 


0.059 


ICU 


0.011 


0.603 


0.016 


0.186 


0.089 


0.007 


0.049 


0.720 


0.223 


MED 


0.001 


0.038 


0.041 


0.008 


0.017 


0.027 


0.084 


0.781 


0.131 


SURG 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


1.000 


0.000 


AMB 


\ 0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 


1.000 


0.000 ) 



M(60) 



This shows that all the states have a fairly low expected visit rate, and most patients would 
transition to an absorbing state with a relatively low number of transitions. 

The final quantities we obtain are the asymptotic state probabilities. These are fairly trivial 
when the process has one or more absorbing states. When this is the case we only need to find 
G(oo). The asymptotic probabilities of beginning in states 1 through 6 and being absorbed in 
states 7 through 9 are: 

ECF HOME DIED 
CCU /0.0575 0.7830 0.1595 \ 
PCCU 0.0615 0.8796 0.0589 
ICU 0.0499 0.7272 0.2229 
MED 0.0840 0.7846 0.1314 
SURG 0.0000 1.0000 0.0000 
AMB \0.0000 1.0000 0.0000/ 

Each row adds to 1, because the probability of reaching an absorbing state at infinity should 
be 1 if the SMP has properly defined distributions for each transition. 

In this example we demonstrated the capabilities of using modern computing power and tech- 
niques with the elegant theory presented in Pyke (1961b). All the computations were done in 



7T 



R and even with transitions that did not have closed-form LTs, the total computation time was 
about 70 seconds (on a computer with a dual core 1.9 GHz processor). The computational time 
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would dramatically decrease if the SMP was parameterized with transition distributions that had 
closed-form LTs. 

We have calculated several quantities that are of use to statistical modelers and are not often 
included in an analysis. Many times only the asymptotic probabilities are reported but these are 
of little use for processes that are not recurrent. 



5 Discussion and conclusions 

In this paper we have presented computational techniques for deriving quantities of interest in 
semi-Markov processes with smooth transition distributions; these techniques can be used in most 
continuous time SMPs. This section describes some performance issues and significant extensions 
which broaden the scope and applicability of SMPs. 



5.1 Performance 

SMPs with many states appear in applications such as queuing problems, where the number of 
states could reach more than 100,000. The methodology we presented works well for SMPs with a 
small number of states, however at some point the number of states will become computationally 
problematic. The computational performance will be affected by a several factors such as the 
sparseness of the q(t) matrix, the number of states, n, and the number of elements in q(s) 
that must be found numerically. In SMPs with a small number of states, we can afford to find 
(I — q(s)) -1 many times, but as n increases it may be beneficial to find it once symbolically using 
programs such as Maple or Mathematica. To calculate SMPs with large n, care must be taken 
in defining states, parameterizing the state transitions and choosing methods for calculation of 
(I-qOO)" 1 . 



5.2 Uncertainty quantification 

Besides solving for quantities such as P i (t) and G i .(t), we are often interested in estimating 
the uncertainty associated with these quantities, based on uncertainties in estimating transition 
probabilities and waiting time distributions. A full treatment is beyond the scope of this paper, but 



we provide a summary here; see (Huzurbazar 2005, Chapt. 5) for further details and examples. 



Writing the waiting time distributions as if the parameters are estimated from 

sampled waiting times, a frequentist uncertainty analysis that accounts for sampling variability 
can be done as follows: draw bootstrap resamples for each transition, re-estimate the parameters 
for each Fij(t;9ij), and solve the process. After repeating this, say, m times, the results can be 
used to construct confidence bands for quantities of interest. 

Alternatively, we can do a Bayesian analysis by assigning prior distributions to the #'s, and 
computing posterior distributions based on waiting time samples. Then solving the process re- 
peatedly using m samples from the posterior of the 0's we obtain m estimates of, say, G.. (t). Using 
these m estimates we can (point-wise in time) find 95% credible intervals for the true G.^t). 
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5.3 Nonparametric and semiparametric analysis 



The methods we have described may be applied to solve SMPs nonparametrically, using sampled 
waiting-time data without assuming parametric distributions. Given a sample x%, X2, ■ ■ ■ , x n of 
waiting times for a given transition, we use the empirical distribution function (EDF) F n (t) and 
the Laplace-Stieltjes transform to define the empirical Laplace transform (ELT) of the data as 



e~ sx dF n {x) 



n 

n ^ 

i=l 



After substitution of fij(s) for fij(s) in the construction of the q matrix, all the algorithms 
described in Section [3] can be used without change to derive, e.g., the empirical transform of a 
first passage distribution; it can be shown (Collins; 2009) that the resulting ELTs are strongly 



consistent estimators of the corresponding exact transforms. Parametric LTs and ELTs may also 
be mixed in constructing q, for a semiparametric analysis. 

Since the EDFs underlying ELTs are not smooth, the EULER method described in Section 



2.2 performs poorly. However, a modified EULER is developed in Collins and Huzurbazar (2009), 



which is both fast and accurate (for large enough samples) for ELT inversion. 



In conclusion, this paper fills a significant gap in the SMP literature by connecting the theory 
with functional formulas that can be used by practitioners in a large variety of applications. 
Our ongoing work is directed towards expanding the applicability of the methods, improving 
performance, broadening the base of computational platforms, and developing software that is 
easier to use. 
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